## Add the estimated vm pop in the drawings to the survey

library(here)
library(tidyverse)
library(sf)
library(sp)

system_info <- Sys.info()
if (system_info["sysname"] == "Darwin") {
  Sys.setenv(PROJ_LIB = "/opt/homebrew/Cellar/proj/9.4.1/share/proj")
}

load(here("Data", "surv_sf_dat.rda"), verbose = TRUE)
load(here("Data", "da_vmpop_avg_lst.rda"), verbose = TRUE)
load(here("Data", "survey.geo.rda"), verbose = TRUE)

## The da_vmpop_avg_lst.rda file comes from polygon_overlap_calculations ---
## maybe later add the names and convert to data frame in that file For now,
## that file runs a long time, even on a machine with 256 GB RAM and 38 cores

names(da_vmpop_avg_lst) <- surv_sf_dat$vcid

## Some people drew maps that did not include any Canadian DA
missing_calcs <- sapply(da_vmpop_avg_lst, function(obj) {
  all(is.na(obj))
})

da_vmpop_avg_lst[missing_calcs]

valid_calcs <- da_vmpop_avg_lst[!missing_calcs]

## Now convert to a data.frame
valid_calcs[1:2]
# $`32407`
#  prop_vmpop_map_avg prop_vmpop_map_dawt
#          0.06970348          0.04725731
#
# $`33504`
#  prop_vmpop_map_avg prop_vmpop_map_dawt
#           0.1858270           0.1607507

calcs_df <- data.frame(do.call("rbind", valid_calcs))
#         prop_vmpop_map_avg prop_vmpop_map_dawt
# 32407         0.0697034754        4.725731e-02
# 33504         0.1858270369        1.607507e-01

calcs_df$vcid <- as.numeric(row.names(calcs_df))

orig_rows <- nrow(survey.geo)
survey2 <- left_join(survey.geo, calcs_df, by = "vcid")
stopifnot(nrow(survey2) == orig_rows)

save(survey2, file = here("Data", "survey2.rda"))
